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NUMERICAL  MODELING  OF  TWO-DIMENSIONAL  WIDTH-AVERAGED  FLOWS 
USING  BOUNDARY-FITTED  COORDINATE  SYSTEMS 


PART  I:  INTRODUCTION 

1.  The  use  of  numerically  generated  boundary-fitted  curvilinear 
coordinate  systems  as  the  basis  for  numerical  solution  of  partial  dif¬ 
ferential  equations  on  arbitrary  regions  is  now  well  established.  A 
comprehensive  survey  of  the  generation  and  use  of  these  coordinate  sys¬ 
tems  has  recently  appeared  (Thompson,  Warsl,  and  Mastin  1983),  and  the 
proceedings  of  a  recent  symposium  devoted  to  this  area  (Thompson  1982a) 
cover  the  basic  techniques  involved. 

2.  Such  coordinate  systems  have  the  property  that  some  coordinate 
line  is  coincident  with  each  segment  of  the  boundary  in  the  physical 
region,  so  that  the  complication  of  boundary  shape  is  effectively 
removed  from  the  problem.  In  the  past  decade,  the  numerical  generation 
of  curvilinear  coordinate  systems  has  provided  the  key  to  the  develop¬ 
ment  of  finite-difference  solutions  of  partial  differential  equations  on 
regions  with  arbitrarily  shaped  boundaries.  Although  much  of  the 
Impetus  for  these  developments  has  come  from  fluid  dynamics,  the  tech¬ 
niques  are  equally  applicable  to  heat  transfer,  electromagnetics,  struc¬ 
tures,  and  all  other  areas  involving  field  solutions. 

3.  With  coordinate  systems  that  make  coordinate  lines  (surfaces 
in  three  dimensions)  coincident  with  the  boundaries,  finite-difference 
codes  can  be  written  which  are  applicable  to  general  configurations 
without  the  need  of  special  procedures  at  the  boundaries.  Even  when  the 
boundaries  are  in  motion,  the  use  of  such  coordinate  systems  allows  all 
computation  to  be  done  on  a  fixed  grid  with  a  uniform  rectangular  mesh 
in  the  transformed  (computational)  plane.  This  greatly  simplifies  the 
coding  with  regard  to  boundary  conditions,  which  can  now  be  represented 
without  need  of  interpolation.  It  is  also  possible  to  distribute  the 
curvilinear  grid  lines  in  the  physical  plane  with  concentration  of  lines 


In  regions  of  high  gradients  while  maintaining  the  square  grid  In  the 
transformed  plane. 

4.  With  such  systems,  the  grid  points  may  be  thought  of  as  a 
finite  set  of  observers  of  the  physical  solution,  stationed  so  as  to  be 
most  effective  in  covering  all  of  the  action  on  the  field.  The  struc¬ 
ture  of  an  Intersecting  net  of  families  of  coordinate  lines  allows  the 
observers  to  be  readily  Identified  In  relation  to  each  other.  This 
results  In  more  simplified  coding  than  would  the  use  of  a  triangular 
structure  or  a  random  distribution  of  points.  The  grid  generation 
system  provides  some  Influence  of  each  observer  on  the  others  so  that 
when  one  moves  to  get  into  a  better  position,  its  neighbors  will  follow 
in  order  to  maintain  smooth  coverage  of  the  field.  The  curvilinear 
coordinate  system  thus  should  cover  the  field,  with  coordinate  lines 
coincident  with  all  boundaries.  The  distribution  of  lines  should  be 
smooth,  with  concentration  In  regions  of  hlgl.  gradient. 

5.  Numerical  solutions  of  partial  differential  equations  are  done 
on  the  curvilinear  coordinate  system  by  first  transforming  all  partial 
derivatives  (or  integrals)  analytically  from  cartesian  to  curvilinear 
coordinates,  so  that  the  latter  become  the  Independent  variables.* 

Normal  and  tangential  derivatives  at  boundaries  are  similarly  trans¬ 
formed.  (These  transformation  relations  are  given  in  Thompson  1982b). 
The  result  is  a  set  of  partial  differential  equations  and  boundary  con¬ 
ditions  in  which  all  derivatives  (or  integrals)  are  taken  with  respect 
to  the  curvilinear  coordinates.  These  equations  may  then  be  expressed 
as  difference  equations  on  a  rectangular  grid  In  the  transformed  plane. 
Thus,  there  is  no  need  for  interpolation,  regardless  of  the  shape  of  the 
boundaries  or  the  distribution  of  the  curvilinear  coordinate  lines  in 
the  field. 

6.  A  finite-difference  solution  is  given  for  the  two-dimensional 
(2D),  time-dependent,  width-averaged  Navier-Stokes  equations,  including 
an  algebraic  turbulence  model,  based  on  such  a  numerically  generated 


A  simple  example  would  be  a  conversion  from  cartesian  to  polar 
coordinates. 


boundary-fitted  coordinate  system.  This  solution  Is  applicable  to  2D 
regions  of  arbitrary  shape,  with  multiple  Inlets  and  outlets,  and  with 
obstacles  In  the  Interior.  The  density  is  taken  to  be  a  function  of  the 
temperature,  and  the  system  of  equations  to  be  satisfied  consists  of  the 
continuity  equation,  the  two  momentum  equations,  and  the  energy 
equation. 

7.  The  solution  provides  for  a  choice  of  central,  upwind,  or  ZIP 
differencing  of  the  convective  terms.  Viscous  terms  and  heat-conduction 
terms  are  represented  by  second-order  central  difference  expressions. 

The  time  derivatives  may  be  represented  as  either  first-  or  second-order 
backward  difference  expressions.  The  finite  volume  formulation  Is  used 
so  that  the  equations  are  fully  conservative.  The  solution  is  implicit 
in  time,  with  all  the  difference  equations  being  solved  iteratively  by 
successive  overrelaxation  in  each  time  step. 

8.  The  input  allows  any  portions  of  the  boundary  (external  or 
obstacles)  to  be  designated  as  inlets,  outlets,  no-slip  surfaces,  or 
slip  surfaces.  Arbitrary  specification  of  velocity  and  temperature  (or 
density)  on  inlets  and/or  outlets  is  allowed.  The  output  is  in  the  form 
of  field  arrays  of  the  velocity  components,  pressure,  and  temperature. 
All  computation  is  done  in  metric  units,  but  the  input  and  output  units 
may  be  specified  otherwise.  The  code  for  this  solution  (WESSEL)  and  the 
code  that  generates  the  boundary-fitted  coordinate  system  (WESCOR)  are 
described  in  detail  in  Thompson  and  Bernard  (1985)  and  Thompson  (1983), 
respectively.  The  coordinate  system  is  generated  from  the  numerical 
solution  of  a  system  of  elliptic  partial  differential  equations  with 
provision  for  controlling  the  spacing  of  the  coordinate  lines  in  the 
field  (Thompson  1982c).  The  transformed  region  is  rectangular,  with  the 
obstacles  and  intrusions  transformed  to  slits  and/or  slabs.  A  complete 
discussion  of  the  relation  between  the  physical  and  transformed  planes 
is  given  in  Thompson  (1983). 

9.  The  mathematical  development  of  this  solution  and  the  finite 
volume  formulation  thereof  on  a  general  boundary-fitted  coordinate  sys¬ 
tem  are  discussed  in  the  following  sections.  This  solution  is  designed 


PART  II:  EQUATIONS  OF  MOTION 


10.  Since  numerical  solutions  are  Inherently  descriptions  of 
physical  properties  Interpreted  either  as  values  at  a  finite  collection 
of  points  or  as  average  values  over  finite  regions.  It  Is  natural  to 
write  the  governing  equations  of  motion  In  the  form  of  Integral  conser¬ 
vation  equations.  In  this  form,  the  equations  express  time  rates  of 
change  In  a  finite  volume  In  terms  of  the  resultant  fluxes  through  the 
boundary  of  the  volume.  Such  equations  make  no  assumptions  regarding 
behavior  Inside  the  volume,  and  naturally  provide  only  average  values  of 
the  solution  over  the  volume. 

11.  It  Is,  In  fact.  In  this  integral  form  that  the  physical  laws 
of  dynamics  and  thermodynamics  must  be  originally  postulated  for  a 
fluid,  since  the  concept  of  density  Is  strictly  definable  only  through 
an  Integral  expression.  Only  If  the  fluid  medium  Is  assumed  to  be  con¬ 
tinuous  can  such  Integral  relations  as  the  Divergence  Theorem  be  used  to 
obtain  partial  differential  equations  by  applying  the  Integral  conserva¬ 
tion  equations  to  an  arbitrarily  small  volume.  In  this  case,  the  Inte¬ 
gral  form  that  represents  the  difference  In  fluxes  on  opposite  sides  of 
the  volume  becomes  a  partial  derivative  at  a  point. 

12.  While  the  Integral  and  differential  forms  are  analytically 
equivalent  for  a  continuum,  the  two  forms  can  have  different  Implica¬ 
tions  In  a  numerical  solution,  which  Is  necessarily  carried  out  on  a 
finite  set  of  points  or  volumes.  This  is  particularly  true  when  the 
solution  uses  points  on  a  non-cartesian  grid  or,  equivalently,  volumes 
that  are  not  cubes.  It  Is  also  true  when  physical  "discontinuities" 
such  as  fronts  and  shocks  occur  In  the  field.  In  each  of  these  cases, 
the  numerical  representation  of  derivatives  may  not  always  be  accurate, 
while  It  may  be  possible  to  represent  fluxes  through  the  volume  sides 
with  sufficient  accuracy. 

13.  For  these  reasons,  the  equations  of  motion  are  given  below, 
first,  in  the  integral  conservation  form.  The  analytically  equivalent 
differential  form  Is  then  given  for  reference  and  some  later  use,  but 


the  developments  of  numerical  procedures  that  follow  are  based  primarily 
on  the  Integral  form.  Throughout  the  presentation  of  the  equations, 
subscripts  Indicate  cartesian  coordinate  directions  In  tensor  notation, 
with  repeated  Indices  Indicating  summation  as  usual.  The  Kronecker 
delta  appears  with  Its  usual  meaning,  as  well.  The  common  symbols  for 
the  physical  variables  are  used.* 


Integral  Form  of  the  Equations 


Continuity  equation 

14.  The  continuity  equation,  expressing  conservation  of  mass  for 
an  arbitrary  moving  volume.  Is 

■A  I  PdV  +  /  /  'p(u,  -  U.)n,dS  -  0  (1) 

Here,  the  first  term  measures  the  rate  of  change  of  mass  in  the  volume, 
while  the  second  term,  being  a  surface  Integral  over  the  entire  closed 
boundary  of  the  volume.  Is  the  resultant  flux  of  mass  out  of  the  volume 
through  Its  boundary.  In  the  second  term,  u^  Is  the  fluid  velocity, 
while  Uj  is  the  local  surface  velocity,  so  that  Uj  -  Uj  Is  the  fluid 
velocity  relative  to  the  surface.  Note  that  (u^  -  simply  the 

dot  product  (u  -  y)*n  and  hence  the  relative  velocity  normal  to  the 
local  Increment  of  bounding  surface,  so  that  the  product 
Is  the  flux  of  mass  through  this  surface  Increment.  If  the  surface  Is 
moving  with  the  fluid,  then  the  relative  velocity,  Uj  -  Uj  ,  vanishes 
and  the  equation  states  simply  that  the  mass  In  the  volume  must  be 
constant . 

Momentum  equation 

15.  Newton's  Second  Law,  that  the  time  rate  of  change  of  momentum 
equals  the  sum  of  the  Impressed  forces,  applies  to  a  particular  mass. 


* 


For  convenience. 


■^endix 


A). 


symbols  are  listed  and  defined  In  the  Notation 
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not  to  an  arbitrary  volume.  Therefore*  the  volume  used  in  this  relation 
must  always  contain  the  same  particles  of  the  fluid,  and  Its  boundary 
surface  must  move  with  the  local  fluid  velocity.  Since  the  fluid 
velocity  relative  to  the  surface  then  vanishes,  the  expression  of  the 
Second  law  Is  as  follows: 


n///. 


Vp(t) 


pu^dV 


K 


s.(t) 


ti^n^dS  + 


///. 


(t) 


Pgj^dV 


(2) 


Here  the  subscript  F  Is  attached  to  the  volume  and  Its  bounding  sur¬ 
face  to  Indicate  that  the  volume  Is  being  defined  to  always  contain  the 
same  portion  of  fluid.  The  left  side  then  is  the  time  rate  of  change  of 
the  component  of  momentum  of  the  particular  mass  of  fluid  within 

the  volume.  The  forces  Impressed  on  this  mass  are  of  two  types,  as 
expressed  by  the  two  terms  on  the  right  side.  The  first  term  represents 
the  forces  acting  on  the  surface  of  the  fluid  volume.  These  forces 
arise  from  pressure  and  viscosity  and  are  expressed  In  terms  of  the 
stress  tensor  (force  per  unit  area)  in  which  the  subscripts  indi¬ 

cate,  respectively,  the  direction  of  the  force  and  the  direction  of  the 
normal  to  the  surface  on  which  It  acts.  The  product  then  is  the 

surface  force  In  x.  direction  acting  on  the  local  surface  Increment 
(i)  ^ 

dS  .  The  second  term  on  the  right  Is  the  body  force,  such  as 
gravity,  acting  on  the  mass  In  the  volume,  the  vector  g^  being  the 
body  force  per  unit  mass.  This  body  force  could  Include  electromagnetic 
forces  as  well  as  gravity,  of  course. 

16.  The  momentum  equation  (Equation  2)  may  be  converted  to  apply 
to  an  arbitrary  moving  volume  through  the  Reynolds  Transport  Theorem, 
which  states  that  for  any  moving  volume  and  any  function  f  : 

T.  -  III,,,  M  ^ 


Applied  to  the  fluid  volume  V_  this  becomes 


fu,n,dS  (3b) 


Now  consider  an  arbitrary  volume  V(t)  that  Instantaneously  is  identi¬ 
cal  to  the  fluid  volume  V  (t)  .  When  the  two  volumes  coincide,  the 

r 

first  terms  on  the  right  sides  of  Equations  3a  and  3b  become  identical 
so  that  we  have  at  that  Instant 


//X  ,,,  -  3?  //(<,,  -  II 


f(Uj  -  U^)njdS  (4) 


Since  any  arbitrary  moving  volume  will  coincide  with  some  fluid  volume 
at  any  instant,  this  last  equation  may  be  used  to  replace  the  left  side 
of  Equation  2,  with  the  substitution  Pu^  ■=  f  .  The  momentum  equation 
may  be  written  for  an  arbitrary  moving  volume  as: 

r.  Ill,,,  "“i-'  ^  ll,„  »“i<“i  - 

■  "  JIl,„ 

Now  the  first  term  on  the  left  is  the  time  rate  of  change  of  momentum  in 
the  volume,  and  the  second  term  on  the  left*  is  the  resultant  flux  of 
momentum  through  the  bounding  surface. 

Energy  equation 

17.  The  energy  equation  expresses  conservation  of  energy  for  an 
arbitrary  moving  volume,  and  can  be  written  directly  as 

5?  ///v(„  4^-2  I  I,,,  K'  "  ^ 

■  II, „  "  IIl,„  ““I*/’  ■  H,o 


"TV-.' ■v’-Tl'v" •A-.’-."'*-''  •7'^' 


Analogous  to  the  other  equations,  the  first  term  on  the  left  represents 
the  time  rate  of  change  of  total  energy  (Internal,  measured  by  tempera¬ 
ture,  and  kinetic)  In  the  volume,  and  the  second  gives  the  resultant 
flux  of  total  energy  through  the  bounding  surface.  The  first  two  terms 
on  the  right  are  the  work  done  by  the  surface  and  body  forces,  respec¬ 
tively,  while  the  last  term  Is  the  resultant  heat  transfer  through  the 
bounding  surface.  In  this  equation,  the  energy  release  by  Internal 
reactions  has  been  omitted.  The  surface  heat  transfer  vector  q^  could 
Include  radiation  as  well  as  conduction. 

Stress  and  heat  conduction 

18.  In  these  equations,  the  stress  tensor  and  the  conduction  heat 
transfer  vector  are  given  by,  respectively. 


Is  the  Kronecker  delta.  Both  of  these  quantities  are 


where 

Inherently  differential,  being  dependent  on  gradients  In  the  fluid,  so 
the  assumption  of  continuity  of  the  medium  Is  necessary  for  their  defi¬ 
nition.  However,  average  values  of  derivatives  can  be  defined  using  the 
Divergence  Theorem  which  states,  for  a  general  vector  A  , 


-/(vj 


dS 


(9) 


If  A  Is  taken  as  k^^^f  ,  where  k^^^  Is  the  unit  vector  In  the 


direction,  then  Equation  9  becomes,  for  any  function  f  , 


//1 87- -/I'M 


dS 


(10) 


which  can  be  used  to  define  an  average  value  of 
terms  of  an  integral  over  the  bounding  surface. 
Summary  of  equations 


over  the  volume  in 


a  i. 


19.  Equations  1,  5,  and  6,  together  with  7  and  8  and  an  equation 
of  state  and  relations  for  the  viscosities  and  conductivity,  constitute 
a  closed  system  of  equations  of  motion  in  integral  form. 

20.  With  the  portion  of  the  stress  tensor  due  to  viscous  forces 
written  as 


these  equations  are  collected  here  as  follows: 


(11) 


Continuity: 


(12) 


Momentum: 


.dS 


III  Pu  dV  +11  Pu.(u.  -  U.)n  c 

=  -  ff  pn.dS  +  ff  a.  n.dS  +  fff  Pg^dV  (13) 


S(t) 


■  'flw  "M(t) 

+  ///  Pu  g  dV-//  q  n  dS  (14) 
-'•'-'V(t)  J  ^  J-'S(t)  ^  ^ 


Energy  equation  In  terms  of  enthalpy 

21.  The  energy  equation  can  also  be  written  In  terns  of  enthalpy, 
rather  than  the  energy,  by  substituting 


h  “ 


e  + 


£ 

P 


to  obtain  for  the  left  side  of  the  energy  equation 


UjuJdV  p(h  +  i 

■  dt  fff  ~fj  ’’*“1 


-  Uj)„jdS 
-  »j)njdS 


Fixed  volume 

22.  For  a  fixed  volume 
inside  the  volume  Integral, 
equation  becomes 


,  Uj  “  0  and  the  time  derivative  passes 
For  example,  the  first  term  of  the  momentum 


3(Pu^) 

3F“ 


dV 


with  similar  changes  In  the  other  equations 


Differential  Form  of  the  Equations 


Basic  equations 

23.  The  differential  form  of  these  equations  is  obtained  by 
applying  the  Divergence  Theorem  (Equation  9)  to  convert  all  surface 
Integrals  to  volume  integrals  and  then  reasoning  that,  since  the  result 
must  apply  for  any  fixed  volume  no  matter  how  small,  the  integrand  must 
vanish.  This  yields  the  following  equations  from  Equations  12,  13, 
and  14  for  the  continuity,  momentum,  and  energy  equations,  respectively. 

Continuity: 


.  3(pu  ) 

^  + _ i_ 

3t  3x. 

3 


=  0 


(15) 


Momentum: 


3t  3x  3x^  9Xj  ^^i 


(16) 


Energy: 


_3 

3t 


[p  (e  +  ^  UiUi)j  +  ^  jp  (e  +  I 


9(pu  )  3(u  0..)  3q. 

- ^  +  pu.g.  -  ^  (17) 

9x.  9x.  ^  9x. 


The  differential  form  of  the  equations  implies  the  assumption  of  con¬ 
tinuity  of  the  medium,  of  course,  since  it  was  necessary  to  apply  the 
Integral  equations  to  an  arbitrarily  small  volume  to  obtain  this  form. 


Jl 


24.  In  the  differential  form,  the  mechanical  energy  contribution 
to  the  energy  equation  (Equation  17)  may  be  removed  by  multiplying  the 
momentum  equation  (Equation  16)  by  u^  and  subtracting  the  result  from 
the  energy  equation,  the  result  being 


9(pe) 


3(peu.) 


ij  9x. 


or  In  terms  of  enthalpy, 


9(ph) 

3t 


3(phu  ) 


3p  .  3e_  .  1 

n  "j  IIT  "ij  35^ 


Thus,  Equations  15,  16,  and  one  of  the  Equations  17-19,  together  with 
equations  of  state,  etc.,  constitute  a  closed  system. 

25.  These  forms  of  the  energy  equation  with  mechanical  energy 
removed  can  be  put  In  Integral  form  by  Integrating  over  the  volume  and 
reversing  the  use  of  the  Divergence  Theorem  (Equation  9)  and  the 
Reynolds  Transport  Theorem  (Equation  3)  with  the  result 


j  j'  j 


"  /i(t)  "'“j  ■ 


dS  (20 


or,  with  the  enthalpy. 


'  ^///»(t,  *Um 

fff  °n  ~  f f 

■^“'■'v(t)  •'-'s(t)  ^  ^ 


Therefore,  in  integral  form,  the  energy  equation  (Equation  14)  can  be 
replaced  by  either  Equation  20  or  21. 

Set  of  eouations  for  nresent  model 


26.  The  set  of  equations  chosen  for  the  present  model  is  composed 
of  the  continuity  equation  (Equation  12),  the  momentum  equation  (Equa¬ 
tion  13),  and  the  energy  equation  (Equation  20),  together  with  the 
stress  and  heat  flux  relations.  Equations  11  and  8,  with  the  bulk 
viscosity  y'  set  to  zero.  These  equations  are  applied  on  a  fixed 
curvilinear  coordinate  system  so  that  the  surface  velocity 
vanishes.  The  time  derivatives  then  appear  as  partial  derivatives 
inside  the  Integrals. 


PART  III:  EQUATIONS  OF  MOTION  IN  GENERAL 
CURVILINEAR  COORDINATE  SYSTEMS 


27.  In  order  to  treat  completely  arbitrary  configurations,  the 
equations  of  motion  must  now  be  transformed  from  the  cartesian  system  of 
the  preceding  section  to  a  general  2D  curvilinear  coordinate  system. 

t 

I 

General  Curvilinear  Coordinates 


Unit  tangents  and  normals 

28.  To  establish  the  terminology,  consider  the  following  general 
element  bounded  by  four  curved  sides,  on  each  of  which  one  of  the  curvi¬ 
linear  coordinates  is  constant: 


line  of 

constant 

(n-line) 


29.  With  the  position  vector  r  =  ix  +  jy  ,  we  have  the  following 
relations.  The  unit  tangent  to  a  line  of  constant  n  is  given  by 


2  2 

with  Y  =  ,  and  the  unit  normal  to  a  line  of  constant  n  is 


n('l)  .  k  X  t^^>  -  (jx,  -  iy,) 
/7 


(2 


Similarly,  for  a  constant-C  line,  the  unit  tangent  and  normal  are 


(O 


dr 

dn 

dr 

dn 

\ 

1 . 

7==  -  ^  ^ 


(2 


n  n 


2  2 

with  a  = 

n  ■'n 


The  normal  is 


^  ^(C)  X  k  >  -p  (-jx^  +  iy^) 


(2 


Area  and  volume 


30.  Then  the  area  of  a  face  on  an  h-line  between  two  5-lines  is 
given  by 


dS 


(n) 


dr 


|k  X  (ix^  +  jy^) I 

-  iy^l  =  +  yl  =  ^ 


(2 


and  the  area  of  a  face  on  a  5-line  between  two  n-lines  is 

dr  I 


dS 


(5) 


(2 


Then  on  an  n-line 


K'- 

i 


^(n)ds(^)  .  -  iy^ 


(28) 


and  on  a  C-llne, 


.  -Jx^  +  ly^ 


Also,  for  any  vector  A  =  iA^  +  jA2  ,  we  have  on  an  n-line. 


A-n^^^dS^"^^  =  -Ajy^  + 


while  on  a  C-line, 

A*n^^^dS^^^  -  A^y^  -  A2X^ 

31.  The  volume  of  a  cell  is  given  by 

I  /dr  dr  \  1 

|k- 

=  |k-[(ix^  +  jy^)  X  (ix^  +  jy^)]| 

=  |k-(kx^y^  -  kx^y^) I 

“  '^C^n  ■  ^ 

where  J  is  the  Jacobian  of  the  transformation. 

Divergence 

32.  Now  by  the  Divergence  Theorem,  we  have  for  any  vector  A 


Application  of  the  Divergence  Theorem  to  the  cell  yields 


-  [A.„'»ds''>] 

Here,  the  following  notation  is  used: 


Thus 

(V.A)J  =  (-A^y^  +  V?^+  ~  Vf^- 


^^^1%  ■  Vn^C+  ■ 


Note  that,  in  the  limit,  this  could  be  written  in  derivative  form 


(V*A)J  =  (-A.y  +  Ax)  + 

~  ~  i  5  C  n 


(Ajy^  -  Vn's 


where  the  subscripts  now  indicate  partial  differentiation. 

33.  This  then  becomes  a  relation  for  the  divergence: 


(33) 


<sa 


V-A  -  i  [(-Aj,^  ^  AjXj)^  +  Ajy„  -  A^x^)  J 


If  the  derivatives  are  expanded,  there  is  a  cancellation  of  cross  derivatives 
of  X  and  y  so  that 


V«A  -  -  l-A,  y.  +  A-  X-  +  A,  y  -  A«  x 
'  '  J  \  ‘y,  ^  2^  c  2^  n 


34.  Both  Equations  33  and  34  are  valid  expressions  for  diver¬ 
gence,  and  the  two  are- equivalent  analytically.  They  are  not  equivalent 
numerically,  however,  since  the  latter  Involves  the  expansion  of  deriva¬ 
tives.  The  form  of  Equation  33  Is  the  geometrically  conservative  form 
(after  multiplication  through  by  the  Jacobian).  Note  that  the  differ¬ 
ence  between  these  two  forms  Is  that  In  the  geometrically  conservative 
form.  Equation  33,  the  "area"  through  which  the  flux  of  A  flows  is 
that  of  the  bounding  faces  through  which  the  flux  occurs.  In  Equation 
34,  however,  the  area  is  evaluated  at  the  cell  center  rather  than  on  the 
boundary.  For  a  sharply  deformed  cell,  such  as  Is  Illustrated  below. 


the  fluxes  through  and  sides  would  both  be  computed  using  the 

area  on  the  dotted  line  in  Equation  34.  With  Equation  33,  however,  the 
actual  areas  of  the  £;+  and  sides  would  be  used.  Thus,  the  fully 

conservative  form  should  be  much  more  tolerant  of  deformed  coordinate 
systems. 


Geometrically  conservative 
derivatives  and  surface  integral 

35.  If  A  is  taken  to  be  if  ,  then  by  Equation  33, 


f 

X 


(35) 


which  gives  the  geometrically  conservative  representation  of  the  partial 
derivative  with  respect  to  x  .  Similarly,  with  A  =  jf  we  have 


(36) 


36.  Recall  also  for  later  use  that  the  right  side  of  Equation  33, 


1 


except  for  the  —  ,  is  the  surface  integral. 


II  II 

JJA.n.dS  =  J  Ja- 


3  3 


ndS  .  (A^y^  -  +  (-A^y^  A  A^x^)^ 


-  '‘P?  +  <*2>n 


(37) 


where 


^1  =  ^1%  “  Vn  ^38a) 

^2  ^  ~^1^C  ^  (38b) 


Here  A^  and  A2  are  the  fluxes  of  A  through  the  C  and  n  faces, 
respectively.  In  fact,  note  that  by  Equations  23  and  25, 


Thus  Aj^  is  the  component  of  A  normal  to  a  C-line  multiplied  by  the 
area  of  the  5-face  (/a).  Similarly,  A  is  the  normal  component  of  A 
on  an  n-line  multiplied  by  the  area  of  the  ii“face. 

37.  Comparing  Equations  33  and  37,  it  is  clear  that  the  equations 
of  motion  can  be  put  in  geometrically  conservative  form  either  by  apply¬ 
ing  the  Integral  form  to  the  cells  of  the  curvilinear  coordinate  system, 
l.e.,  using  Equation  37,  or  by  simply  representing  the  derivatives  in 
the  differential  form  by  the  geometrically  conservative  expression, 
l.e.,  using  Equation  33. 

Geometrically  noncon- 
servative  derivatives 

38.  In  geometrically  nonconservative  form  we  have,  with  the 
derivatives  expanded  in  Equations  35  and  36, 

f  -  4  (y  fr  -  yrf  )  (39) 

X  J  •'n  K  ^  n 


f 

y 


I 

J 


V?' 


(40) 


39.  The  second  derivatives  are,  by  repeated  application  of  these 
relations. 


f 

XX 


u 


- 


)-7^( 


DDYY  f  +  DDYX  f 


(41) 


40.  Finally t  the  Laplaclan  ist  from  Equations  41  and  42 

■  7  ^ ’O  -  7  ^  fx) 


with 


DDY  5 


DDX  =  -  26Xj^  +  yx^^ 


6  5  XjX^  +  y^y^ 


and  a  and  y  defined  earlier: 


a 


Y 


Transformed  Equations  of  Motion 


Continuity  equation 

41.  Using  Equation  37  we  have 


UjnjdS 


(pu), 


+ 


where  by  Equation  38 


Here  u  and  v  are  the  velocities  normal  to  the  C  and  n  lines, 
respectively,  multiplied  by  the  area  of  the  faces  through  which  the  flow 
occurs  (/a  and  /y  »  respectively).  Then,  since 


the  continuity  equation  (Equation  12)  becomes 


p^J  +  (pu)^  +  (pv)^  *  0 


(46) 


As  noted  above,  this  form  could  also  have  been  obtained  by  using  the 
geometrically  conservative  expressions  for  the  derivatives,  i.e.. 
Equations  35  and  36,  in  the  differential  form  of  the  continuity  equation 
(Equation  15).  The  same  applies  for  the  geometrically  conservative 
forms  of  the  momentum  and  energy  equations  developed  below. 

Momentum  equation 

42.  In  the  momentum  equation  (Equation  13)  we  have,  by  Equation  37 


// 


dS  -  (pu^u)^  +  (pu^v)^ 


lj..jdS  -  (o^j)^  + 


where,  by  Equation  38 


The  pressure  term  Is  a  bit  different.  We  have 


jjpadS  - 

.  [p5<«>ds<«)]j.  -  L<'>dS<5>]j. 

■  [p<i»5  -  *  [^<4%  *  i%)] 


by  Equations  28  and  29.  Then 


Jjppds  .  i[(py„)^  -  (pyj)^  •  J  [-(P%)5  *  (p-j) 


Finally,  the  gravity  term  becomes 


lljpg^dV  -  pg^J 


The  entire  momentum  equation  then  Is 


(pu)^J  +  (puu  -  “  ®12  ”  “ 


(pv)  J  +  (pvu  -  a,, 


-  px^)_  +  (pvv  -  a,,  +  px  )  -  pg,J 


Energy  equation 

43.  For  the  energy  equation  (Equation  20)  we  have,  by  Equation  37 


// 

J  Jpeu.n 


j..jdS  “  (peu)^  +  (pev)^ 


Ik: 


dS  -  (qj)j  +  (,2)„ 


with,  by  Equation  38 


qj  =  -  q2% 


^2  “  ‘‘2*5  ■ 


The  energy  equation  then  becomes. 


dUj  3u 


(pe)^J  +  (peu  +  qp^  +  (pev  +  q^)^ 


(50a) 


(50b) 


‘'^ij  ^  -  P  ^  °  ° 

J  J, 


(51) 


Stress  and  heat  conduction  terms 


\  •  3  [<“%>{  - 

(52a) 

"y  ■  3 

]  (52b) 

'x  ■  J  ['''^>{  - 

(52c) 

''y  ■  3 

]  (52d) 

"k  ■  3  [‘*Ve  -  <">'e>n] 

(52e) 

^  ■  3  [-<%>£  * 

]  (52f) 

With  these  relations,  we 

then  calculate  from  Equation  11,  with  y’  -  0  , 

^11  =  3  "“x  -  1 

(53a) 

<^22  "  3 

(53b) 

°12  “  <^21  '  ‘‘S 

(53c) 

and  from  Equation  8, 

‘’l  ■  -  -^^x 

(54a) 

<^2  ’  -  '^■^y 

(54b) 

Also,  for  the  last  term 

In  the  energy  equation  (Equation  51)  we  have 

31 


Wall  heat  transfer 


45.  When  a  cell  side  coincides  with  a  wall,  the  wall  heat 
transfer  supplies  the  value  of  q  on  that  side  as  follows: 


For  an  n-llne  wall 
normal  as  given  In  Equation  23, 


we  have,  using  the 


S„ll  ■  %■  ^ 


±1  ~ 


where  the  +  applies  to  the  lower  wall  and  the  -  to  the  upper.  T 


^2  '  -  'Iwall 


A 


A 


A 


46.  Similarly,  for  a  ^-llne  wall 


(C) 


4„<«) 

/ 


we  have 


q^  -  ±  /a 


S/all 


with  the  +  corresponding  to  the  left  wall  and  the 


to  the  right 


Wall  temperature 

47.  For  an  n-llne  wall,  =  ±  S/all  Equation  56.  Then 

using  Equations  50b,  54a-b,  and  52e-f  we  have 

*  S,all  ■  ’2  ■  Vs  -  ‘■iS'S 

•  <V«  ■  ■'x>'s) 

- r  |[-<N>S  ‘Vn]^  •  [‘%>5  ■  <Vn]>'s} 

or, 

(Tx^)^x^  +  (Ty^)^y^  -  ^  ‘‘wall 

The  wall  temperature  can  be  evaluated  from  this  equation  using  one-sided 
differences  for  the  n  derivatives. 

48.  For  a  5-line  wall, 

(Tx^)5X^  +  (Ty^)5y^  -  (Tx5)^x^  +  (Ty5)^y^  +  ^  q^^^^  (59) 

Here  one-sided  5  derivatives  are  used. 

49.  With  the  derivatives  expanded,  the  nonconservative  forms  of 
these  conditions  are,  from  Equations  58  and  59,  respectively, 

rT„  .  6l^  ;  i  (60) 


for  an  l-line,  wall  and 


'*wall 


(61) 


aT^  -  6T  +  - 
C  n  ic 


for  a  5-line  wall. 

Slip  boundary  conditions 

50.  On  a  boundary  with  free  slip,  the  normal  velocity  component 
and  the  vorticity  must  vanish.  These  two  conditions  serve  to  determine 
the  two  velocity  components  on  such  a  boundary  as  follows. 

51.  The  vorticity  is  given  by 

0)  =  V  -  u  =0  (62) 

X  y 

or  in  curvilinear  coordinates,  using  Equations  52a“C, 

(vy  ).  -  (vy,)  +  (ux  )  -  (ux  )  “  0  (63) 

nC  Cn  nC  Cn 

52.  On  an  n~line  boundary,  the  vanishing  of  the  normal  velocity  is 
expressed  from  Equation  45b  as 


vx^  -  uy^  =  0  (64) 

and  in  the  above  expression  of  vanishing  vorticity,  we  have,  in  discrete 
form,  with  the  upper  sign  applying  to  the  lower  wall  w  ,  etc. 

where  (F^,  F2,  F^)  are  (1,  -1,  0)  and  -2,  for  first  and  second 

order,  respectively.  Then 


Thus 


0 

-  ''3p''’'e>w±2  " 


Then,  with 


r'"*  =  *  ['''y„'£  *  '“%>£]-  '■2 [<’'’'£>«!!  *  ‘"■‘e'wsi] 

■  ^3[‘''''{>wS2  *  <"’‘£>.±2] 


ve  have 


u 


.  Jk  £<”> 

3^1 


V 


p(n) 


33.  Similarly  on  a  5-line  boundary,  we  have,  by  Equation  45a 


so  that  with 


we  have 


-  ’'3  ['''^n’vty  *  <“n>„i2] 


-  Jjl 

aF 


(67a) 


aF 


Model  Set 


(67b) 


54.  The  set  of  equations  in  general  curvilinear  coordinates  used 
in  the  present  model  are  the  continuity  equation  (Equation  46) ,  the 
momentum  equation  (Equation  48),  the  energy  equation  (Equation  51), 
together  with  the  expressions  given  by  Equations  52-55  and  the  appli¬ 
cable  boundary  conditions  given  by  Equations  56  and  57  for  the  heat 
transfer,  Equations  60  and  61  for  the  wall  temperature,  and  Equations  65 
and  67  for  the  wall  slip  velocity.  Additional  relations  and  some  modi¬ 
fications  to  these  equations  are  presented  in  Part  IV. 
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PART  IV;  MODIFICATIONS  TO  THE  EQUATIONS  OF  MOTION 


Removal  of  Hydrostatic  Pressure 


55. 

pressure 

equation 


One  modification  Is  made  by  subtracting  out  the  hydrostatic 

p„  which  satisfies,  by  Equation  13  with  zero  velocity,  the 
n 


-  iL^n^dS  +  IK 


where  p  Is  the  Initial  density.  The  momentum  equation  (Equation  13) 
s 

then  can  be  written,  for  a  fixed  volume,  as 


.dS 


///^  ffvfi 

*  -  //<P  -  *  fll^”  ■  P,**! 


dV  (68 


with  the  hydrostatic  pressure  calculated  from 


=  Po  j 


(69 


Here  p^  is  the  pressure  at  some  reference  point 

56.  This  modification  is  reflected  in  the  momentum  equation 
(Equation  48)  by  the  replacement  of  the  pressure  with  the  difference 
between  the  true  pressure  and  the  hydrostatic  pressure,  and  replacement 
of  the  density  in  the  gravity  term  by  the  difference  from  the  Initial 
value.  This  Is  valid  only  if  the  fluid  is  initially  in  a  condition  of 
static  equilibrium. 


57.  As  a  stabilizing  measure,  the  products  of  the  discrete  resid¬ 
ual  of  the  continuity  equation  (Equation  12)  and  the  velocity  and  energy 
are  added  to  the  right  sides  of  the  momentum  and  energy  equations. 
Equations  48  and  51,  respectively,  as  corrective  source  terms.  These 
terms  are,  respectively,  u^D  and  eD  ,  where  D  is  the  residual 
evaluated  from  Equation  46: 


°  =  ///ft  *  //“j"j‘ 


=  p^J  +  (pu)^  +  (Pv)^ 


58.  It  may  be  noted  that  the  analytical  effect  of  this  corrective 
use  of  the  continuity  equation  in  the  momentum  and  energy  equations 
would  be  to  reduce  the  equations  in  the  nonconservative  form  if  the 
equations  were  coverted  to  differential  form  and  terms  were  cancelled 
where  possible.  However,  the  effect  is  not  the  same  in  the  difference 
expressions  since  the  use  of  different  points  in  the  various  expressions 
prevents  such  cancellations  on  the  discrete  grid.  The  numerical  effect 
of  these  corrective  terms  is  to  resist  the  tendency  of  large  outflow 
from  a  cell  to  catastrophically  deplete  the  cell  contents,  and  similarly 
to  lessen  the  effects  of  large  inflow. 

Lateral  Averaging 


59.  Variable  width  is  accounted  for  in  the  present  model  by  the 
procedure  of  lateral  averaging,  which  has  seen  much  use  in  hydrodynamic 
problems  to  introduce  some  three-dimensional  effects  into  two- 
dimensional  computational  models.  With  this  procedure,  the  equations  of 
motion  are  integrated  laterally,  and  the  dependent  variables  are 


replaced  by  lateral  averages  thereof.  Thus,  with  the  z-direction  being 
the  lateral  direction,  we  have  the  definition  of  the  lateral  average  of 
a  quantity  f  given  as 


f (x,y) 


1 

B(x,y) 


/■ 


z,(x,y) 


f (x,y,z)dz 


Zj(x,y) 


(7i; 


where  B  Is  the  width,  l.e.,  B(x,y)  =  Z2(x,y)  -  Zj(x,y)  .  Then,  with 
f  a  perturbation  from  this  average  value,  we  have 


f(x,y,z)  =  f(x,y)  +  f*(x,y,z)  (Jl\ 

In  all  that  follows,  such  perturbations  will  be  assumed  to  be  small  In 
comparison  with  the  average  value.  Thus  |f'|  <<  |fl  . 

Continuity  equation 

60.  Applying  this  procedure  to  the  continuity  equation 
(Equation  15),  we  have 


Neglecting  the  quadratic  prime  terms,  we  then  have  for  the  laterally 
averaged  continuity  equation 


'  It  *  al  97 

The  effect  of  the  lateral  averaging  on  tne  continuity  equation  In  geo¬ 
metrically  conservative  form  (Equation  46)  Is  thus  to  multiply  the 
density  by  the  width  B  In  all  terms. 

Momentum  equation 

61.  Applying  the  same  procedure  to  the  x-momentum  equation  from 
Equation  16,  we  have,  with  subscripts  Indicating  components,  not 
derivatives. 


similar  terms  occur  for  the  other  convective  terms.  The  brackets  cancel 
as  in  the  continuity  equation.  The  convective  terms  then  become 


3(BpU^)  3(BpUv)  r  r  rin  n 

— ^ +  prime  terms 

Neglecting  the  p'  terms,  we  obtain  for  the  remaining  prime  terms: 


62.  Also,  for  the  pressure  terms. 


Now  consider  the  pressures  impressed  on  the  volume: 


If  p  and  p-  are  written  p  -  p  -  ^  (p  +  p..,)  ,  then 


66.  In  the  present  model,  the  effect  of  the  correlation  terms  is 
Incorporated  Into  the  averaged  stress  terms  through  turbulence  models  as 
discussed  in  a  later  section.  The  stresses  on  the  lateral  boundaries 
are  neglected.  The  hydrostatic  pressure  Is  removed  from  the  above  equa¬ 
tions  by  simply  defining  the  pressure  in  the  above  equations  to  be  the 
difference  between  the  true  pressure  and  the  hydrostatic  pressure  and 
replacing  the  density  in  the  gravitational  terms  with  the  difference 
from  the  reference  value. 

67.  In  the  geometrically  conservative  forms  of  the  momentum  equa¬ 
tion  (Equation  A8) ,  the  effect  of  the  lateral  averaging  is  to  multiply 
the  first  terms,  the  pressure  terms,  and  the  gravity  terms  by  the  width 
B  ,  and  to  Insert  B  in  the  ^  and  n  derivatives  for  the  convective 
and  stress  terms. 

Energy  equation 

68.  In  a  similar  manner,  the  lateral  averaged  energy  equation  is, 
using  Equation  18, 

B  +  a(Bpue)  ^  3(Bpve)  ^  -  j 3(Bu)  ^  3 (By) [ 

3t  3x  3y  P  L  3*  J 

f-  3(Bu)  -  (  3(Bv)  3(Bu)\  -  3  (Bv)~| 

j_°xx  3x  ®xy  \  3x  9y  /  °yy  9y  J 


3(Bq  )  3(Bq  ) 

+  — =  0 


3x 


3y 


(75) 


The  correlation  terms  between  velocity  and  energy  have  been  neglected, 
along  with  the  stress  work  and  heat  transfer  on  the  lateral  boundaries. 

69.  The  effect  of  the  lateral  averaging  on  the  energy  equation  in 
geometrically  conservative  form  (Equation  51)  is  to  multiply  the  first 
term  by  the  width  B  and  to  insert  a  B  inside  the  ^  and  n  deriva¬ 
tives  and  inside  the  derivatives  in  the  last  term. 


Effect  on  hydrostatic 
pressure  and  mass  residual 

70.  Since  the  width  B  Is  multiplied  by  the  pressure  and  gravity 
terms  In  the  laterally  averaged  momentum  equation  (Equation  74),  the 
lateral  averaging  has  no  effect  on  the  hydrostatic  pressure  In 
Equation  69. 

71.  The  laterally  averaged  mass  residual  Is  obtained  by  inference 
from  Equations  70  and  73.  Dropping  the  bars  Indicating  average,  and 
using  the  transformed  derivatives,  this  quantity  becomes 

D  =  Bp^J  *  (Bpu)^  +  (Bpv)^  (76) 

Incompressible  Flow 

72.  Of  particular  Interest  in  the  present  work  are  situations  In 
which  the  density  variations  are  negligible  in  all  terms  except  the 
gravity  term  pg  .  Attention  Is  further  restricted  to  cases  in  which 
density  may  vary  with  temperature  but  not  with  pressure.  (That  is,  the 
temperature  gradients  are  large  enough  to  affect  the  density,  but  the 
pressure  gradients  are  not.)  Thus,  from  a  practical  standpoint,  the 
fluid  is  mechanically  incompressible  and  the  laterally  averaged  con¬ 
tinuity  equation  reduces  to 

(Bu)  +  (Bv)  =  (Bu)c  (Bv)  =  0 

X  y  5  n 


which  represents  the  well-known  Boussincsq  approximation  for  variable- 
density  flow.  Note  that  this  removes  the  pressure  term  from  the 
laterally  averaged  energy  equation  (Equation  75). 


Final  Equations 


73.  The  final  set  of  equations  to  be  solved  In  the  present  model 
then  consists  of  the  momentum  equation,  the  energy  equation,  and  the 
continuity  equation  (discussed  later)  together  with  the  relations  for 
the  stress  and  heat  conduction  and  the  boundary  conditions.  This  set  is 
summarized  as  follows: 

Momentum  equation 

74.  From  Equation  48  with  the  modifications  Indicated  above  we 

have 


Bp 


11  J  +  B(p  uu  -  a,  , )  +  1 

o  t  [_  o  Jf  L 


B(PoUV  -  o^2> 


-In 


+  B(P^y^  -  -  B(p  -  Pg)gjJ  -  uD  =  0  (77a) 


and 


BPoVtJ  +  [b(p^vS  -  ^21)]^  +  [bCp^vv  -  022) 


1 


+  B(-P^x^  +  P^x^)  -  B(p  -  Pg)g2J  -  vD  =  0  (77b) 

where  P  e  p  -  is  the  difference  from  the  hydrostatic  pressure,  p^ 
Is  the  initial  density,  p  is  the  present  density,  and  p^  is  a  con¬ 
stant  reference  density.  Nonconservative  expressions  have  been  used  for 
the  pressure  terms,  for  reasons  that  will  be  explained  later. 


Energy  equation 

75.  From  Equation  51  with  the  modifications  we  have 


BPoCtJ  +  ‘’rj  '^2'J 


3(Bu  ) 

-  J  a,  .  — ^  -  eD  =  0  (78) 

Ij  3*j 


Continuity  equa¬ 
tion  and  mass  residual 


76.  Subject  to  the  Boussinesq  approximation,  the  continuity  equa¬ 


tion  is 


(Bu)^  + 


and  the  mass  residual  is  now  given  by 


■  [<“=){  •  <®'>n] 


Boundary  conditions 


77.  The  boundary  conditions  are  as  follows: 

£.  Wall  heat  transfer:  Equations  56-57. 

Wall  temperature:  Equations  60-61. 

£.  Slip  wall  velocity:  Equations  65  and  67. 

78.  The  following  choices  of  boundary  conditions  are  incorporated: 

£.  Thermal. 

(1)  Specified  boundary  temperature,  with  the  boundary 
heat  transfer  calculated  from  Equations  56  or  57. 

(2)  Specified  heat  transfer,  with  boundary  temperature 
calculated  from  Equations  60  or  61. 

(3)  Extrapolated  boundary  temperature  from  interior, 
l.e.,  equal  to  the  interior  value  adjacent  to  the 
boundary. 

£.  Velocity. 

(1)  Specified  boundary  velocity. 


■•V'  .> 


(2)  Slip,  with  the  boundary  velocity  calculated  from 
Equations  65  and  66. 

(3)  Extrapolated  boundary  velocity  from  interior. 

Solid  walls  are  treated  with  thermal  and  velocity  boundary  condition  (1) 
or  (2).  Inlets  have  thermal  and  velocity  boundary  condition  (1)  or  (3). 
Outlets  have  thermal  and  velocity  boundary  condition  (1)  or  (3).  Free 
surfaces  are  approximately  simulated  with  thermal  and  velocity  boundary 
condition  (3). 

79.  The  normal  velocities  on  all  inlets  and  outlets  must  be  such 

that  the  total  outflow  exactly  balances  the  total  inflow.  The  mass  flow 

rate  through  an  n-llne  boundary  is,  by  Equation  45  with  cognizance  of 

the  width  B  ,  simply  while  that  through  a  C-llne  boundary  is 

p^Bu  .  Outflow  through  a  lower  (upper)  n-llne  boundary  is  then  -(+) 

p  Bv  ,  and  Inflow  is  +(-)  p  Bv  .  Similarly,  outflow  through  a  left 
o  o 

(right)  5-llne  boundary  is  -(+)  »  while  inflow  is  +(-)  P^Bu  •  The 

total  flow  rates  are  calculated  by  summing  over  all  points  on  the  appro¬ 
priate  boundary  segments. 

Auxiliary  relations 

80.  In  the  above  equations,  the  following  auxiliary  relations  are 


£.  u  and  V  from  Equation  45. 

Jb.  from  Equation  49. 

£.  q^  from  Equation  50. 

from  Equation  53. 

£.  q^  from  Equation  54. 

Velocity  and  temperature  derivatives  for  stress  and  heat 
conduction  from  Equation  52. 

£.  The  relation  Equation  55,  with  u  and  v  multiplied  by 

B  . 


h^.  Pjj  from  Equation  69. 

D  from  Equation  81. 

In  addition,  the  conservative  derivative  expressions  are  used  for  the 
Xj-  derivatives  in  the  energy  equation.  The  metric  components  are 
repeated  here  for  completeness: 


PART  V:  DISCRETE  REPRESENTATION 


81.  The  present  model  is  a  finite-difference  solution  based  on 
the  integral  form  of  the  equations  of  motion,  i.e.,  with  the  equations 
in  geometrically  conservative  form.  The  solution  is  implemented  on  a 
boundary-fitted  coordinate  system  which  allows  boundaries  of  arbitrary 
shape  to  be  treated. 

82.  The  boundary-fitted  coordinate  system  is  generated  numer¬ 
ically  by  solving  a  system  of  elliptic  partial  differential  equations  as 
discussed  in  general  in  Thompson  (1982c),  and  in  detail  for  this  present 
application  in  Thompson  (1983).  The  latter  reference  also  discusses  the 
operation  and  input  of  the  coordinate  code  used.  The  essential  feature 
of  this  curvilinear  coordinate  system  is  that  some  coordinate  line  seg¬ 
ment  is  coincident  with  each  segment  of  the  physical  boundary.  Including 
Interior  obstacles.  This  allows  all  difference  expressions  to  be  taken 
along  coordinate  lines  regardless  of  the  boundary  shape. 

83.  V/lth  such  a  coordinate  system,  all  computation  can  be  done  on 
a  square  grid  on  the  rectangular  transformed  region  regardless  of  the 
shape  of  the  physical  region.  This  allows  the  physical  configuration  to 
be  changed  via  the  input  without  modification  to  the  code. 

Solution  Grid 


84.  The  coordinate  system  is  generated  with  twice  as  many  points 
in  each  direction  as  is  intended  for  use  in  the  flow  solution.  Thus, 
the  grid  system  for  the  flow  solution  is  formed  by  taking  every  other 
coordinate  line.  The  physical  variables  are  all  defined  at  the  grid 
points  on  this  grid.  Values  needed  between  grid  points  are  determined 
by  simple  linear  averages.  However,  since  the  coordinate  values  are 
available  at  points  between  these  grid  points,  coordinate  derivatives  at 
points  between  the  grid  points,  as  well  as  at  the  grid  points,  may  be 
calculated  without  averaging  the  coordinate  values.  This  is  Important 
in  achieving  a  conservative  difference  form,  since  averaging  the 


coordinate  values  can  cause  different  sets  of  coordinate  values  to  be 
used  in  equivalent  representations  of  certain  derivatives  with  a  conse¬ 
quent  Introduction  of  spurious  effective  gradients  in  uniform  functions. 


Difference  Representation 

85.  The  points  surrounding  a  point  of  calculation,  l.e.,  point  C 
in  Figure  1,  are  Identified  In  relation  to  the  compass  directions.  As 
noted  above,  the  coordinate  derivatives  at  any  of  the  25  points  on  this 
figure  can  be  calculated  directly  from  the  coordinate  values  at  adjacent 
points,  e.g.. 


(x^)c  =  Xg  - 
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since  the  coordinate  system  is  generated  using  twice  as  many  points  in 
each  direction  as  will  be  used  in  the  flow  solution. 

86.  The  flow  solution,  however.  Is  represented  only  at  grid 
points,  so  that  averages  must  be  used  at  other  points.  Thus 
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Figure  1.  Coordinate  scheme  for  finite-difference  representation 

With  these  values  one  can  calculate  the  derivatives  of  the  flow  variable 
at  all  points,  as  Illustrated  for  the  coordinates  above. 

87.  The  flow  variables  are  stored  in  2D  arrays,  so  that  the  fol¬ 
lowing  correspondencies  apply: 


-  F(I.  J) 


fm  ■  fd.  J  ’  » 


^NENE  -  F(I  •  1.  J  ♦  1).  etc. 


88.  The  volume  of  integration  used  is  that  indicated  by  dotted 
lines  in  the  figure,  and  the  Jacobian  J  thus  represents  the  area  of 
the  four-sided  area  in  physical  space  corresponding  to  this  square  in 
the  transformed  plane.  The  fluxes  must  therefore  be  represented  on  the 
sides  of  this  square  in  the  transformed  plane,  meaning  that  the  ^  and 
n  derivatives  in  the  momentum,  energy,  and  mass  residual  equations 
(Equations  77,  78,  and  80)  are  represented  as 


Difference  Equations 


89.  The  finite-difference  representations  for  the  various  terms 
in  the  equations  are  discussed  below.  All  types  of  terms  are  covered, 
but  the  full  equations  are  not  given  here  because  the  code  (Thompson  and 
Bernard  1985)  is  written  with  all  terms  Identified  therein. 

Time  derivatives 

90.  All  time  derivatives  are  represented  by  first-  or  second- 
order  backward  difference  expressions: 


where  the  superscript  indicates  the  time  level  and  the  coefficients 
(Cj,  c^,  c^)  are  (1,  -1,  0)  for  the  first  order  and  -2,  for  second 
order.  Second  order  is  used  for  all  but  the  first  step,  where  the  first 
order  must  be  used  to  start  the  solution. 

Convective  terms 

91.  The  convective  terms  are  of  the  form  (Bp  fu).  or  (Bp  fv)  , 

o  c  on 

where  f  is  1  for  continuity,  u  and  v  for  momentum,  and  e  for 
energy.  These  terms  can  be  represented  by  central  or  upwind  differ¬ 
encing,  the  choice  being  designated  in  the  input.  Central  differences 
are  second  order  while  upwind  is  first  order.  Although  the  conservative 

form  Bp  f  J  +  (Bp  fu)_  +  (Bp  fv)  is  standard,  it  is  also  possible  to 
O  t  O  5  O  f) 

choose  the  nonconservative  form  Bp  Jf  +  Bp  uf)  +  (Bp  vf)  obtained 

o  t  o  ^  on 

by  using  the  continuity  equation  in  the  momentum  and  energy  equations, 
or  the  ZIP  form  (Zalesak  1981),  via  an  input  parameter. 

92.  The  conservative  central  difference  expressions  are  of  the 

form 

while  the  nonconservative  central  forms  are 

-  V  <«») 
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The  conservative  upwind  forms  are 


[(BpoUf) -  (Bp^a)„f^  (84 

[(BpoVf) Jc  =  (BPo^>N^T  -  ^®Po^^S^B 

where 


The  nonconservative  upwind  forms  are 


<»C<,='j)c  '  «'>0=>C»R  - 

with  the  definitions  of  the  subscripts  R  ,  L  ,  T  ,  and  B  as  above 
except  that  all  the  decisions  are  based  on  the  values  u  and  v_  . 
93.  The  ZIP  form  is 


where  here  the  term  (Bp^uf)j^  represents  the  sum  of  four  terms,  in  each 
of  which  one  of  the  four  variables  Is  evaluated  at  point  E  and  the 
other  three  are  evaluated  at  point  C  .  The  other  terms  in  Equation  86 
have  analogous  meanings,  the  subscripts  L  ,  T  ,  and  B  being  asso¬ 
ciated  with  points  W  ,  N  ,  and  S  ,  respectively,  in  the  same  manner. 

94.  Regardless  of  the  mode  selected  for  the  convective  terms,  the 
central  form  is  always  used  for  the  evaluation  of  mass  residual 
(Equation  80)  in  the  momentum  and  energy  equations. 

Stress  and  heat  conduction  terms 

95.  The  stress  terms  are  evaluated  using  second-order  central 

differences  as  follows.  The  term  Equation  77a  is  given  by 


(Boii)^  =  (BojPj.  -  (Ba^j)„ 


Now,  by  Equation  49a, 


(3;j)e  -  (cr„y^  - 


and,  by  Equation  53, 
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Finally,  by  Equation  52a, 


and  analogous  expressions  for  the  other  terns  Involved,  lliis  pattern  is 
followed  for  all  the  stress  terms  In  the  momentum  equation. 

96.  In  the  energy  equation  (Equation  78),  the  dissipation  term 
3(Bu^) 

a..  — - -  Is  evaluated  in  a  similar  fashion.  Thus,  for  this  term, 

IJ  3Xj 

which  only  occurs  at  point  C  ,  we  have 
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Also,  for  this  term,  the  derivatives  I — - -  are  evaluated  as  is 

1  C 

done  for  (u  )„  above.  '  •'  ' 

X  C 

97.  The  heat  conduction  terms  In  the  energy  equation  (Equa¬ 
tion  78}  are  evaluated  analogously  to  the  stress  terms  in  the  momentum 
equations,  using  second-order  central  differences.  Thus 


with,  by  Equation  50a, 


<’l>E  ■  <‘'l^  -  Vh’e 


and,  by  Equation  54a, 
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with 

«x>E  ■  [<%>5  -  <^c>n]E 

and 
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etc . 

Pressure  terms 

98.  The  pressure  terms  in  the  momentum  equations  are  of  the  fol' 
lowing  forms,  which  are  evaluated  using  first-order  backward 
differences : 


<'’c  -  W 
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The  reason  for  the  one-sided  differencing  is  explained  below. 


Chorin's  method  for 
calculating  pressure 


99.  For  incompressible  flow,  the  pressure  gradient  serves  one 
purpose:  it  prevents  loss  or  gain  of  mass  within  a  fluid  volume.  Thus, 
the  pressure  and  the  continuity  equation  are  directly  related.  In  curvi 
linear  coordinates,  the  desired  condition  for  conservation  of  mass  is 
given  by  Equation  79,  which  is  the  same  thing  as  saying 

D  =  0 

with  the  mass  residual  D  given  by  Equation  80.  The  question,  of 
course,  is  how  to  find  a  pressure  distribution  that  maintains  this 
condition. 

100.  Consider  the  divergence  of  the  width-averaged  momentum  equa¬ 
tion,  obtained  by  summing  the  x-derlvatlve  of  Equation  48a  and  the 
y-derivatlve  of  Equation  48b: 

“t  *  *  4  ^ 

Note  that  all  terms  not  involving  the  pressure  or  the  mass  residual  have 

been  placed  anonymously  in  the  right-hand  side  (RHS) .  The  expressions 

for  P  and  P  are  the  nonconservative  forms  given  by  Equations  39 
X  y 

and  40.  Employing  a  backward  time-differencing  scheme  for  ,  Equa¬ 
tion  89a  can  be  rewritten  as 

H  *  a?  <%'’x  -  4 

with  D  representing  the  current  mass  residual,  RHS  containing  the 
values  of  D  from  previous  time  steps,  and  c^  being  the  leading 
coefficient  in  Equation  81. 

101.  Now  suppose  that  the  left-hand  side  (LHS)  of  Equation  89b  is 
calculated  by  a  convergent  Iteration  scheme  such  that 


where  the  superscript  (m)  indicates  the  iteration  count.  If  the 
scheme  is  truly  convergent,  the  condition  D  =  0  can  be  achieved  by 
setting 


=  0 

in  each  successive  calculation  of  the  left  side  of  Equation  89c.  This 
in  itself  defines  an  iteration  scheme,  namely 


Equation  89d  is  a  Poisson  equation  for  pressure  that  must  be  solved  at 
each  iteration  for  the  entire  flow  field.  What  remains  is  to  select  the 
difference  expressions  that  will  represent  the  quantities  ,  P^  , 
and  D  . 

102.  In  the  numerical  solution  of  the  momentum  equations  (Equa¬ 
tions  77a  and  77b),  one  can  choose  upwind,  central,  or  ZIP  differencing 
for  the  advective  terms,  with  central  differencing  for  the  viscous 
terms.  Experience  has  shown,  however,  that  central  differencing  for  the 
pressure  terms  in  the  momentum  equation,  and  for  the  mass  residual  in 
the  continuity  equation,  promotes  numerical  oscillation  on  regular 
grids.  The  way  to  avoid  the  latter  problem  is  to  use  a  staggered  grid 
(pressure  and  velocity  calculated  at  alternate  grid  points)  or,  for  the 
present  work,  to  use  one-sided  differencing  for  the  pressure  gradient 
and  for  the  continuity  equation.  While  this  does  not  reduce  the 


accuracy  of  the  momentum  calculation.  It  does  reduce  the  conservation  of 
mass  to  a  first-order  approximation. 

103.  Assuming  the  coordinate  derivatives  constant  to  first  order, 
but  retaining  B  as  a  variable,  the  pressure  terms  In  Equation  89b 
reduce  to 


al  '  1  a!  «"■{>  -  7  a! 

4  («XjPy  -  ByjP^)  =.  1  ^  (BP^)  -  f  4  (BPj) 


The  one-sided  difference  approximations  for  the  mass  residual  in  Equa¬ 
tion  89d,  and  for  the  pressure  derivatives  In  Equations  82a,  82b,  and 
89d,  are  as  follows 


These  approximations  reduce  Equation  89d  to 


(0T)<")  -  ^  *  ,^B^)p; 
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where  OT  represents  pressure  terms  other  than  the  central  pressure 

104.  The  final  step  in  defining  the  iteration  scheme  for  pressure 
is  to  equate  the  OT  terms  on  each  side  of  Equation  89e,  such  that 


C  C 


2At(acBE  -  ^  Y^Bj^) 


[(Bi)EE  -  (Bi)E+  (Bv)jjj,  -  (Bv)^ 


Numerical  experimentation  has  shown  that,  for  nonsquare  grids,  this 
scheme  converges  only  with  underrelaxation,  l.e.. 


(m)  _  (m-1) _ _ 

C  C  2At(aj,BE  - 


[(bG)..  -  (bG)„+  (bG)„„  - 


Equations  89f  and  89g  constitute  a  variation  of  Chorln's  (1967)  method 
for  calculating  pressure.  The  scheme  is  applied  on  boundary  points  and 
field  points  alike. 

Boundary  conditions 

105.  In  all  expressions  for  temperature  and  velocity  boundary 
conditions,  second-order  central  differences  are  used  for  derivatives 
along  the  boundary,  and  first-order  one-sided  differences  are  used  for 
derivatives  off  the  boundary.  Thus,  for  example,  on  an  n-line  boundary 
with  the  field  to  the  north,  we  have 


(90a: 


(90b: 


Corners 

106.  Convex  corners  in  the  transformed  plane  are  treated  by 
averaging  the  results  of  separate  applications  of  the  temperature  and 
velocity  boundary  conditions  on  the  two  faces  of  the  corner. 

107.  Concave  corners  are  treated  simply  by  averaging  the  extra¬ 
polated  values  along  the  two  faces  of  the  corner,  using  first-  or 
second-order  extrapolation  as  selected  on  the  input. 


Calculation  Procedure 

108.  The  momentum  and  energy  equations  are  solved  simultaneously 
at  each  time  step  by  successive  overrelaxation  (SOR)  iteration.  In  the 
difference  representation,  all  coefficients  of  values  at  the  central 
point  C  are  factored  (i.e.  grouped)  together  to  speed  up  the  iterative 
convergence.  In  the  energy  equation,  this  requires  a  short  inner 


iteration  since  the  energy  and  temperature  at  the  central  point  both 
appear  in  the  equation.  The  SOR  iteration  involves  first  calculating 
new  values  from  the  difference  equations  and  then  taking  the  new  itera¬ 
tive  values  of  velocity  and  temperature  as  weighted  averages  of  these 
values  and  the  values  from  the  previous  iteration,  the  weight  assigned 
to  the  newly  calculated  value  being  the  acceleration  parameter  u  .  The 
field  is  swept  repetitively,  calculating  new  velocities,  temperature, 
and  pressure  at  each  point  in  succession.  Before  each  field  sweep,  the 
boundary  values  are  updated.  These  iterative  sweeps  are  continued  until 
convergence  occurs  or  until  a  prescribed  number  of  Iterations  has  been 
completed.  In  the  latter  case,  the  solution  may  be  directed  to  stop  and 
store  the  partially  converged  results  or  to  continue  to  the  next  time 
step.  Output  is  provided  in  the  form  of  plotted  velocity  vectors  and 
contours  of  density,  temperature,  or  pressure,  as  well  as  printed  field 
values . 

Acceleration  Parameters 

109.  The  acceleration  parameters  for  the  SOR  iteration  are  calcu¬ 
lated  locally  for  the  momentum  and  energy  equations  (see  Thompson, 
Thames,  and  Mastln  1977).  With  these  equations  represented  by  the 
general  form 


where  D  (not  to  be  confused  with  the  mass  residual  here)  contains  all 

_ J _ J _ _  A _ _ J  _ - _ _ _ _ _ _ 1  _ ^  J  ^  J  1  __ 


0) 


2 


(928) 


if  4A.  <  B.  and  4A.  <  B 


2  2 
4a2  -bJ 


|C  -  2(A^  +  A^)! 


4A^ 

-B? 

2 

2 

I  +  ly  ic  -  2(Aj  +  A^)!  +  1 


In  all  other  cases,  the  underrelaxation  represented  by  the  second  of 
these  situations  is  used.  Here  1  and  J  are  the  field  diiq^nsious. 

110.  In  the  momentum  equations  (Equations  77a  and  77b)  we  have 
neglecting  derivatives  of  the  width  B  and  defining  F  =  4/3  , 
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C  “  Jp 


where  G  is  1  for  first-order  time  derivatives  and  y  for  second-order 
By  similar  developments,  we  have  for  the  energy  equation  (Equation  78) 


“  A2  -  3  (a  ♦  y) 
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Turbulence  Model 

111.  The  cumulative  effect  of  turbulence  and  the  correlation 
terms  from  the  lateral  averaging  is  modeled  by  using  horizontal  and 
vertical  eddy  viscosity  and  diffusivlty  coefficients.  This  is  done  be 
replacing  the  viscosity  p  by  u  pD,,  (no  summation)  in  the  stress 


,  and  the  conductivity 


by 


In  the  heat 


elements 

flux  c]  •  *^11  **12 

viscosities,  respectively.  Three  variations  of  Dj^j^  and  D22  are  pro- 


Thus ,  D  and  D 


K  +  P  -  Cj 

are  the  horizontal  and  vertical  eddy 


vided  for  input  selection:  (a)  uniform  values;  (b)  the  model  of  Edlnger 
and  Buchak  (1979),  as  given  in  Johnson  (1981);  and  (c)  the  model  of  Kent 
and  Pritchard  (1959),  as  given  in  Brandsma  and  Dlvoky  (1976). 

112.  The  model  of  Edlnger  and  Buchak  is  Implemented  in  the  pre¬ 
sent  model  as 


D 


11 


0.16 


11 


n  +  lOR, 


>  0 


^22  ~  ^ 


V. 


2At 


R^  j<  0 


Cj  -  1.33  Djj 


+  3.33  rJ 


3/: 


R^  >  0 


^2  “  ^22 


R  <  0 


The  model  of  Kent  and  Pritchard  is  implemented  as 
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113.  In  both  of  these,  the  Richardson  number  is  given  by 


and  d  is  the  local  depth,  d,.  is  the  total  depth,  and  h  is  the 

u  I  .  I  - 1  I 


height  above  the  bottom,  i.e.,  h  *  d«  -  d  .  Also  if 


is  less 


than  0.7  — ,  the  former  is  reset  to  the  latter  value. 

114.  The  model  also  includes  an  optional  artificial  eddy 
viscosity  and  conductivity  in  the  form 


"ii 


C  =  c  -  -liLl 

1  2  pR 


where  h  is  the  mass  residual  (Equation  8U)  calculated  with  centr.al 
dif ferences . 


PART  VI:  RESULTS 


115.  The  discrete  representation  for  the  equations  of  motion, 
with  the  solution  thereof,  has  been  implemented  in  the  computer  code 
WESSEL  (Thompson  and  Bernard  1985).  Boundary-fitted  curvilinear  coordi¬ 
nate  systems  are  generated  numerically  by  a  separate  code,  WESCOR 
(Thompson  1983),  allowing  boundaries  of  arbitrary  shape  to  be  used. 

116.  The  WESSEL  code  was  run  for  three  comparison  cases  with 
distinct  geometries: 

£.  Selective  withdrawal  from  a  rectangular  channel  with 
linear  stratification. 

_b.  Cold  Inflow  into  a  flume  with  nonuniform  width  and  depth. 

£.  Flow  over  a  crested  weir,  with  and  without  stratification 
Willie  this  range  of  geometries  and  flow  conditions  Is  not  exhaustive,  it 
does  afford  a  sampling  of  typical  applications  for  the  present  model. 

117.  In  the  calculations  discussed  below,  conservative  upwind 
differencing  was  used  for  the  convective  terms,  with  second-order  back¬ 
ward  dlf ferencirig  for  the  time  derivatives.  Variable  acceleration 
parameters  were  used  for  the  momentum  and  energy  equations,  with  the 
boundary-temperature  acceleration  parameter  set  at  unity.  The  accelera¬ 
tion  parameter  for  the  Chorln  pressure  calculation  was  fixed  at 

w  =  0.88 

A  value  closer  to  unity  might  have  Improved  convergence  in  some  cases, 
but  the  chosen  value  proved  satisfactory  for  all  cases  and  geometries 
considered.  Calculations  were  made  without  using  artificial  viscosity 
and  conductivity,  and  the  turbulence  models  were  inactive  unless  other¬ 
wise  noted.  As  previously  stated,  the  coordinate  systems  are  generated 
v.'ith  twice  as  many  points  as  the  flow  fields.  The  Initial  velocity 
field  was  always  specified  such  that  the  initial  flow  was  horizontal, 
with  uniform  flow  rate. 

118.  Results  are  presented  primarily  in  the  form  of  velocity 
vectors  and  temperature  (or  density)  contours  at  selected  time  steps. 


All  runs  were  made  on  the  CRAY-1  computer  at  the  Boeing  Company 
Bellevue,  Wash. 


Selective  Withdrawal  from  a  Rectangular 
Channel  with  Linear  Stratification 

119.  A  pair  of  classic  papers  by  Pao  and  Kao  (1974)  and  by  Kao, 
Pao,  and  Wei  (1974)  quantified  the  mechanisms  associated  with  the  estab 
llshment  of  selective  withdrawal  from  reservoirs.  The  first  comparison 
calculation  In  the  present  work  duplicates  the  physical  conditions  for 
one  of  the  experiments  reported  In  Kao,  Pao,  and  Wei  (1974). 

120.  Consider  a  rectangular  channel  of  uniform  width  and  depth, 
with  linear  density  stratification,  such  that 

F  -  0.014 


where  F  la  the  denslmetrlc  Fourde  number  given  by 


F  - 


with 

q  *  flow  rate* 

N  -  the  Valsala  frequency 
d  ■  channel  depth 


N  - 


The  specified  value  of  F  was  achieved  herein  by  setting 


*  In  this  case  (linear  stratification),  the  flow  Is  vertically  sym¬ 
metric.  The  values  of  Q  and  d  specified  here  represent  half  of 
the  total  flow  rate  and  vertical  depth. 


g  “  32.2  ft/sec 


3  3 

=  1  gm/cm  =  1.94  slug/ft 

^  =  0.0071  —82—  =  0.0138  slug/ft^ 

<iy  c.  3 

ft-cm 

w  =  width  =  1  ft 
Q  *  0.00684  ft^/sec 

The  experimental  channel  length  was  30.5  ft,  but  for  the  calculation 
this  was  reduced  to  23  ft,  which  was  long  enough  for  the  comparison  of 
results  near  the  outlet.  The  coordinate  grid  used  was  101  x  41  mesh, 
with  a  51  X  21  flow  field  mesh  (Figure  2),  representing  the  top  half  of 
a  vertically  symmetric  flow  field.  Figure  2  shows  only  the  last  4  ft  or 
so  of  the  channel.  The  longitudinal  grid  spacing  Is  uniform  (0.1  ft) 
for  the  last  3  ft,  and  the  vertical  spacing  Is  uniform  (0.05  ft)  every¬ 
where  In  the  flow  field.  Detailed  results  were  desired  only  for  the 
last  3  ft,  so  exponential  spacing  was  used  for  the  20  ft  upstream. 

121.  The  reference  or  characteristic  time  Is  the  inverse  Valsala 
frequency. 


n“^  -  2.09  sec 

and  the  associated  nondimens Iona 1  time  is 

t*  -  Nt 


f  factors  foi 
)  units  Is  pi 


0.209  sec  (At*  =  0.1).  Further  reduction  of  the  time  step  produced  no 
significant  change  in  the  results.  The  number  of  Iterations  was  fixed 
at  15  per  time  step. 

123.  The  Inflow  boundary  condition  (1  1)  was  extrapolated 

velocity  and  density  with  flow  balance  employed.  The  upper  boundary 
(J  =  21)  was  idealized  as  a  no-slip  wall  with  fixed  density.  The  sym¬ 
metry  boundary  (J  =  1)  was  treated  as  a  slip  wall  with  fixed  density. 

The  downstream  boundary  (i  =  51)  was  modeled  as  a  no-slip  wall  with 
extrapolated  density,  except  for  the  outlet  (1  *  51 ,  j  =  1  to  4)  where 
the  density  was  extrapolated  and  the  velocity  was  held  fixed.  The  Ini¬ 
tial  condition  was  a  uniform  horizontal  flow  field  (u  =  0.00684  ft/sec). 

124.  Figure  3  shows  comparisons  of  calculated  and  measured  veloc¬ 
ities  versus  time  along  the  symmetry  line  (j  =  1).  Stations  1,  2,  and  3 
are,  respectively,  1,  2,  and  3  ft  upstream  from  the  outlet.  Certain 
facts  are  readily  obvious: 

£.  The  velocities  at  early  times  are  correctly  predicted. 

_b.  The  steady-state  velocities  are  underpredicted. 

£.  Accuracy  Improves  with  distance  from  the  outlet. 

The  smaller  time  step  (At*  =  0.1)  provides  better 
resolution  of  the  transient  behavior,  but  the  larger 
(At*  »  1.0)  achieves  essentially  the  same  accelerations 
(somewhat  delayed)  and  steady-state  velocities. 

The  deterioration  of  the  symmetry-line  velocity  predictions  with  time, 
near  the  outlet,  can  be  understood  In  part  after  an  examination  of  the 
predicted  and  observed  steady-state  velocity  profiles.  Figure  4  shows 
the  calculated  and  observed  profiles  at  stations  1-3.  The  sharpest  pro¬ 
file  occurs  nearest  the  outlet,  where  the  highest  acceleration  occurs  on 
the  symmetry  line.  The  numerical  calculation  underpredicts  the  maximum 
velocity  primarily  because  the  symmetry  line  (j  ■>  1)  Is  idealized  as  a 
slip  wall,  which  means  the  horizontal  velocity  is  extrapolated  from  the 
next  line  up  (j  =  2) ,  so  the  code  does  not  resolve  the  large  velocity 
gradient  between  j  °  1  and  j  -  2  .  Hence,  as  the  slope  of  the  ob¬ 
served  profile  becomes  gentler  upstream  of  the  outlet,  the  agreement 
between  calculation  and  experiment  Improves. 


STATION  1 


STATION  3 


Figure  4.  Nondlmensional  steady-state  velocity 
profiles  (t*  -  60)  for  selective  withdrawal 
from  rectangular  channel 


125.  Figure  5  presents  vector  plots  of  the  developing  flow  field 
between  station  3  and  the  outlet.  This  shows  clearly  the  growing  region 
of  circulation  at  t*  “  10  and  t*  *  20  ,  and  the  plot  for  t*  ■  60 
represents  the  steady  state.  The  results  In  Figure  5  were  obtained  with 
At*  -0.1 

Cold  Inflow  Into  a  Flume  with  Nonunlform  Width  and  Depth 

126.  The  experimental  facility  represented  In  this  comparison  Is 
the  Generalized  Reservoir  Hydrodynamics  (GRH)  flume  at  WES,  described  In 
Johnson  (1981).  The  flume  Is  80  ft  long  with  a  3-  by  3-ft  cross  section 
at  the  downstream  end.  The  cross  section  of  the  upstream  end  Is  1  by 

1  ft,  and  the  width  expands  uniformly  with  length  over  the  first  20  ft, 
to  a  cross  section  1  ft  deep  and  3  ft  wide.  Thereafter,  the  width 
remains  constant  while  the  depth  grows  uniformly  with  length,  to  3  ft  at 
the  downstream  end.  Plan  and  side  views  of  the  GRH  flume  appear  In 
Figure  6. 

127.  The  water  In  the  flume  was  at  rest  and  homogeneous  at  the 
beginning  of  the  test,  with  the  temperature  being  70.6*  F.  Cold  water 
was  Input  1.5  ft  from  the  upstream  end  at  a  temperature  of  62.0*  F.  A 
baffle  restricted  the  cold  water  to  enter  through  the  lower  0.5  ft  of 
the  cross  section.  The  Inflow  rate  as  0.022  cfs,  with  the  outflow  rate 
at  the  downstream  end  being  the  same.  The  outflow  was  removed  from  a 
l-ln.-dlam  port  located  0.5  ft  above  the  bottom  of  the  flume  and  1.5  ft 
from  either  side.  The  flow  was  observed  to  be  relatively  smooth  and  was 
probably  laminar,  or  at  least  only  In  the  transition  range  to  turbulent 
flow.  Numerical  results  from  several  codes  are  compared  with  experi¬ 
mental  results  In  Johnson  (1981). 

128.  For  simulation  with  the  WRSSEL  code,  a  33  25  coordinate 

grid  was  used,  with  a  17  x  13  flow  field  grid  (Figure  7).  The  inlet  was 
specified  as  the  lower  0.5  f t  (1  ■  1 ,  j  ■  1  to  7)  at  the  upstream  end. 
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Figure  6.  Plan  and  side  views  of  GRH  flume 

and  the  outlet  was  idealized  as  a  single  grid  space  (i  ■  17,  j  ■  3) 

0.5  ft  from  the  bottom  at  the  downstream  end.  Velocities  were  held 
fixed  at  the  inlet  and  outlet,  with  the  initial  (horizontal)  velocity 
field  specified  such  that  the  flow  rate  was  0.022  cfs  everywhere  in  the 
flume.  The  temperature  was  fixed  at  the  inlet  and  extrapolated  at  all 
other  boundaries.  The  free  surface  was  modeled  as  a  slip  boundary,  and 
the  no-sllp  condition  was  employed  on  the  walls. 

129.  Calculations  were  made  for  (a)  laminar  flow,  (b)  the  Kent 
turbulence  model,  and  (c)  the  Edlnger  turbulence  model.  In  each  case, 
the  time  step  was  set  at  5  sec  with  a  limit  of  50  Iterations  per  time 
step,  and  the  code  was  run  for  500  time  steps. 

130.  Figure  8  shows  a  comparison  of  the  VESSEL  results  with 
experimental  data  and  with  predictions  of  the  LARM  code  of  Edlnger  and 
Buchak  (1979).  The  underflow  front  position  is  correctly  predicted  by 
the  laminar  calculation  and  by  the  calculation  using  the  Kent  turbulence 
model.  The  total  drop  in  outlet  temperature  is  underpredicted,  possibly 


Figure  7.  Boundary-fitted  coordinate  system  for  GRH  flume 

due  to  the  coarse  grid  and  rather  crude  treatment  of  the  outlet  Itself 
Furthermore,  while  the  rate  of  temperature  drop  Is  Initially  correct, 
the  onset  Is  probably  delayed  by  the  diffusive  effects  of  upwind  dif¬ 
ferencing.  Evidence  for  this  can  be  seen  In  the  calculation  with  the 
Edlnger  turbulence  model,  which  Is  far  more  diffusive  than  the  Kent 


TIME,  MIN 


Figure  8.  Comparison  of  results  for  cold 
underflow  in  GRH  flume 

model.  Here  the  large  eddy  viscosity  of  the  Edlnger  model  severely 
retards  the  motion  of  the  underflow  front  as  well  as  the  onset  of 
temperature  drop  at  the  outlet. 

131.  Flow  field  development  with  time  (velocity  vectors  and 
temperature  contours)  for  the  three  flume  calculations  can  be  seen  in 


Figures  9-14.*  The  laminar  case  (Figures  9  and  10)  Is  almost  Indis¬ 
tinguishable  from  that  with  the  Kent  turbulence  model  (Figures  11 
and  12).  The  cold  water  flows  along  the  bottom  to  the  outlet,  causing 
recirculation  In  the  warm  water,  with  a  cooling  effect  that  eventually 
works  Its  way  upstream.  The  Edlnger  turbulence  model  (Figures  13 
and  14)  produces  results  that  are  quite  different  from  the  first  two 
cases.  Instead  of  flowing  just  along  the  bottom,  the  cold  water  dif¬ 
fuses,  forming  a  thicker  and  slower-moving  underflow  front.  This  delays 
the  temperature  drop  at  the  outlet  but  causes  more  cooling  In  the  flume 
as  a  whole. 

Flow  over  a  Crested  Weir 

132.  As  a  final  example,  consider  a  rectangular  channel  of  uni¬ 
form  width  (1  ft)  with  a  crested  weir  (Figure  15).  Flow  tests  were  con¬ 
ducted  for  this  configuration  to  verify  a  finite-element  model.**  In 
one  test,  the  flow  (0.074  cfs)  was  initially  stratified  upstream  of  the 
weir,  but  unstratified  downstream.  In  another  test,  the  flow  (0.7  cfs) 
was  unstratlfled  throughout  the  channel.  Figure  16  shows  the  density 
profile  for  the  stratified  case,  4  ft  upstream  of  the  weir  crest. 

133.  WESSEL  calculations  were  made  for  (a)  constant  density  with 
total  discharge  Q  =  0.7  cfs  ,  and  (b)  initial  density  stratification  on 
both  sides  of  the  weir  with  Q  =  0.074  cfs  .  In  the  latter  case,  the 
initial  density  profile  was  that  given  in  Figure  16.  The  boundary- 
fitted  coordinate  system  (Figure  17)  consisted  of  a  101  x  35  grid,  with 
a  51  X  18  flow  field  grid.  The  no-slip  condition  was  employed  on  the 
channel  walls,  and  the  free  surface  was  idealized  as  a  slip  boundary. 
Density  was  extrapolated  at  all  boundaries.  The  outlet  was  0.6  ft  high 

*  The  velocity  vectors  along  the  bottom  are  actually  parallel  to  the 
bottom.  They  do  not  appear  so  in  the  vector  plots  because  the  flume 
Is  not  drawn  to  scale.  This  also  steepens  the  apparent  horizontal 
temperature  gradients  In  the  contour  plots. 

**  Personal  Communication,  1983,  J.  L.  Grace,  Jr.,  Hydraulics  Labora¬ 
tory,  US  Army  Engineer  Waterways  Experiment  Station,  Vicksburg,  Miss. 
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Figure  12.  Temperature  contours  for  calculation  with  Kent 
turbulence  model  in  GRH  flume  (contour  interval  ■  0.5°  F) 
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Figure  13.  Velocity  vectors  for  calculation  with  Edinger 
turbulence  model  in  GRH  flume  (Sheet  1  of  3) 
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Figure  14.  Temperature  contours  for  calculation  with  Edinger 
turbulence  model  in  GRH  flume  (contour  interval  ■  0.5*  F) 

(Sheet  1  of  3) 
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Figure  17.  Boundary-fitted  coordinate  system  for 
channel  with  crested  weir 

(j  =  1  to  7)  at  the  downstream  boundary  (1  =  51),  and  the  outflow  veloc¬ 
ities  were  held  fixed  and  uniform  in  both  cases.  The  Inlet  was  taken  to 
be  the  entire  upstream  vertical  boundary  (1.7  ft  high,  1  -  1,  j  ®  1 
to  18).  The  Inflow  velocities  were  held  fixed  and  uniform  for  the 
unstratlfled  calculation.  In  the  stratified  case,  however,  the  inflow 


velocities  were  extrapolated,  using  the  flow  balance  option,  to  accommo¬ 
date  possible  Inflow  variations  caused  by  density  changes  upstream  of 
the  weir.  The  time  step  for  both  cases  was  set  at  At  =  0.1  sec  ,  and 

the  number  of  iterations  was  limited  to  25  per  time  step  for  the  first 

50  steps,  and  15  per  time  step  thereafter. 

134.  Figures  18  and  19  present  comparisons  of  the  measured*  and 

calculated  velocity  profiles,  at  the  upstream  toe  of  the  weir,  for  the 

unstratified  and  stratified  cases,  respectively.  The  calculated 
unstratlfled  profile  shows  a  stronger  velocity  gradient  than  the  data 
near  the  bottom,  because  the  inlet  condition  (uniform  velocity)  was 
specified  too  close  to  the  weir  (4.125  ft  from  weir  crest)  for  a 
parabolic  profile  to  develop.  The  stratified  calculation  (Figure  19) 
exhibits  better  agreement,  but  the  backflow  in  the  dense  layer  near  the 
bottom  was  not  observed  experimentally.  Again,  this  is  probably  due  to 

(a)  the  proximity  of  the  computational  inflow  boundary  to  the  weir,  and 

(b)  the  extrapolated  inflow-velocity  boundary  condition.  Improved 
results  might  be  obtained  by  adding  a  longer  upstream  section  to  the 
grid,  with  fixed  inlet  velocity  and  density, 

135.  Figure  20  shows  a  comparison  of  calculated  and  observed 
velocity  profiles  3.5  ft  downstream  of  the  weir  crest,  for  unstratlfled 
flow.  The  data  are  taken  from  the  measurements  for  Q  =  0.074  cfs  , 
since  no  downstream  results  are  available  for  Q  =  0.7  cfs  .  The 
calculated  results  represent  a  unit  discharge  of  0.074  cfs,  but  they 
exhibit  the  same  general  behavior  as  the  data  for  0.7  cfs. 

136.  Direct  comparisons  of  velocity  vector  plots  are  shown  for 
the  stratified  and  unstratified  calculations  in  Figures  21-23.  These 
illustrations  show  clearly  the  dramatic  contrast  between  the  flow  pat¬ 
terns  for  the  two  cases.  Figure  24  shows  time-lapse  pictures  of  the 
developing  circulation  downstream  of  the  weir  crest  for  the  unstratified 


*  Personal  Communication,  J.  L.  Grace,  Jr.,  op.  cit. 
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Figure  20.  Velocity  profile  for  unstratified  flow  3.5  ft 
downstream  of  weir  crest 
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Figure  21.  Comparison  of  velocity  vectors  for  stratified  and 
unstratlfled  flow  In  channel  with  crested  weir 


106 


T4  -W 


60  Seconds 


Unstratified 


0.7  cfs 


60  Seconds  Stratified  0.074  cfs 


Figure  22.  Comparison  of  velocity  vectors  for  stratified 
and  unstratlfled  flow  upstream  of  weir  crest 
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Figure  24.  (Concluded) 


flow.  Figures  25  end  26  present  the  same  sequence  of  events,  on  both 
sides  of  the  weir  for  the  stratified  case.  As  pointed  out  above,  the 
backflow  in  Figure  25  is  apparently  caused  by  the  specified  inflow- 
velocity  boundary  condition.  Density  contours  for  the  stratified  flow 
appear  in  Figure  27. 


10  Seconds 


Figure  25.  Velocity  vectors  for  stratified  flow 
upstream  of  weir  crest  (Continued) 
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PART  VII:  CONCLUSION 


137.  The  VESSEL  code  Is  applicable  to  two-dimensional,  time- 
dependent,  width-averaged,  variable-density  flow  in  arbitrary  regions 
with  multiple  Inlets,  Internal  obstacles,  and  boundary  intrusions. 
Changes  in  the  physical  configuration  can  be  made  easily  via  input.  At 
present,  however,  the  code  is  not  capable  of  simulating  true  free  sur¬ 
face  behavior. 

138.  VESSEL  solves  the  Navier-Stokes  equations  without  assuming 
the  pressure  to  be  purely  hydrostatic,  and  this  feature  sets  it  apart 
from  most  hydraulic  codes  now  In  use.  Thus,  VESSEL  provides  a  needed 
addition  to  the  mathematical  modeling  capability  of  the  Corps  of  Engi¬ 
neers,  offering  a  means  of  simulating  fully  convective  flow  easily  and 
economically.  Potential  applications  Include  selective  withdrawal, 
spillway  approach  flow,  outlet  works  design,  and  pump  station  analysis. 

139.  The  accuracy  of  the  time-dependent  predictions  can  be 
expected  to  Improve  with  further  experience  in  the  use  of  the  code. 

Since  VESSEL  is  written  to  be  very  general  in  application,  its  operation 
is  slowed  somewhat  by  this  generality.  As  it  now  stands,  about  1  sec  of 
CRAY  computer  time  (scalar  mode)  is  required  to  execute  15  Iterative 
sweeps  on  a  50  X  20  grid.  A  faster  code  of  more  limited  scope  could  be 
obtained  by  deleting  a  number  of  contingency  checks  in  the  present  ver¬ 
sion.  Further  development  and  evaluation  is  planned  along  these  lines. 

A  new  version  of  the  code,  using  stream  function  and  vorticity,  has 
already  been  developed  and  will  be  documented  in  the  future. 

140.  The  one-sided  differencing  scheme  for  the  continuity  equa¬ 
tion  arose  from  a  discussion  at  VES.*  This  scheme  allows  the  use  of  a 
regular  grid,**  while  retaining  the  point-to-point  velocity  and  pressure 
coupling  of  a  staggered  grid.t  A  more  extensive  investigation  of  the 
method  is  anticipated  in  future  work. 

*  Personal  Communication,  1983,  R.  S.  Bernard,  J.  F.  Thompson,  and 
P.  J.  Roache,  US  Army  Engineer  Vaterwsys  Experiment  Station, 
Vicksburg,  Miss. 

**  Pressure  and  velocity  calculated  at  the  same  grid  points, 
t  Pressure  and  velocity  calculated  at  alternate  grid  points. 
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APPENDIX  A;  NOTATION 


Lateral  width 
Mass  residual 

Internal  (thermal)  energy  per  unit  mass 
Denslmetrlc  Froude  number 
Arbitrary  function 
Gravitational  acceleration 
Components  of  gravity  vector 
Enthalpy  per  unit  mass 

Unit  vectors  In  x-  and  y-dlrectlons »  respectively 

Valsala  frequency 

Components  of  surface-normal  vector 
Surface-normal  vector 
Superscript  Indicating  time-step  number 
Pressure 

Hydrostatic  pressure 
Reference  pressure 
Flow  rate 

Components  of  heat  flux  vector 

Position  vector 

Surface  of  moving  volume  V(t) 

Temperature 

Time 


Nondlmcnslonal  time 


Components  of  surface  velocity 

Components  of  fluid  velocity 

Fluid  velocity  vector 

Nondlmenslonal  velocity 

X-  and  y-components  of  u  ,  respectively 

Fluxes  of  u  through  surfaces  of  constant  ^  and  constant 
respectively 

Components  of  cartesian  position  vector 
Cartesian  coordinates 


2^2 


Kronlker  delta 
Gradient  operator 
Laplaclan  operator 
Incremental  operator 
Dynamic  viscosity 
Bulk  viscosity 
Thermal  conductivity 
Density 

Initial  density 
Reference  density 
Components  of  stress  tensor 
Components  of  shear  stress  tensor 
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